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Abstract 

A  useful  representation  of  signal  data,  besides  being  a  complete  and  stable 
transformation  of  the  information,  should  make  explicit  useful  features  in  the 
data.  In  computer  vision,  the  one-parameter  family  of  images  obtained  from 
the  Laplacian  of  a  Gaussian-filtered  version  of  the  image,  parameterized  by  the 
width  of  the  Gaussian,  has  proven  to  be  a  useful  data  structure  for  the  extrac- 
tion of  feature  data.  In  particular,  the  zero-crossings  of  this  so-called  scale- 
space  data  arc  associated  with  edges,  and  were  proposed  by  Marr  and  others  as 
the  basis  of  a  representation  of  the  image  data.  The  question  arises  as  to  wheth- 
er the  representation  is  complete  and  stable.  We  survey  some  of  the  results  and 
studies  related  to  these  questions,  and  survey  several  papers  that  attempt  recon- 
structions based  on  this  or  related  representations.  We  then  formulate  a  new 
method  for  the  reconstruction  from  zero-crossings  in  scale-space,  based  on 
minimizing  equation  error,  and  present  results  showing  that  the  reconstruction 
is  possible,  but  can  be  unstable.  We  further  show  that  the  method  applies  when 
gradient  data  along  the  zero-crossings  is  included  in  the  representation,  and 
demonstrate  that  the  reconstruction  is  then  stable. 


1.  Notions 

1.1.  Representations 

In  the  fields  of  signal  analysis  and  image  processing,  the  first  stage  of  a  complete  pattern 
recognition  system  typically  applies  some  numerical  process  to  the  digital  data.  In  image  pro- 
cessing, for  example,  it  is  generally  considered  useful  to  extract  edges,  comers,  and  textured 
regions  in  the  image.  In  other  words,  the  features  and  salient  symbolic  infonnation  about  signal 
and  image  data  are  dynamically  associated  with  groups  or  regions  of  the  data,  and  typically 
depend  more  explicitly  on  primitive  features  such  as  discontinuities,  extrema,  and  local  statistical 
properties  of  the  spatially-sampled  data.  It  seems  reasonable,  therefore,  to  transform  the  initial 
data  to  intermediate  representations  to  make  the  primitive  features  more  accessible  to  the  algo- 
rithms for  signal  analysis.  The  collection  of  all  intermediate  representations,  which  might 
include  binary  edge  images,  texture  measures,  and  other  feature  detectors,  comprise  a  representa- 
tion that  replaces  the  original  signal  for  analysis  purposes.  This  is  the  central  idea  for  example,  in 
vision  processing,  of  Marr's  "Primal  Sketch"[l]  or  Tenenbaum  and  Barrow's  "Intrinsic 
Images. "[2] 

Since  all  analysis  is  done  on  the  intermediate  representation,  the  representation  should  carry 
all  the  information  necessary  for  the  inteipretation  of  the  data.  Of  course,  the  main  idea  is  that 
the  relevant  information  should  be  more  explicit  than  the  original  samples,  and  that  data  redun- 
dancies should  be  eliminated.  The  representation  may  constitute  a  data  compressed  version  of 
the  original  sampling,  but  might  as  easily  contain  more  bulk  data,  in  the  attempt  to  represent  dif- 
ferent features.  It  should  be  possible  to  reconstruct  a  version  of  the  original  signal  from  the  inter- 
mediate representation  since  the  representation  should  contain  all  the  essential  information. 
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While  the  reconstruction  need  not  be  mathematically  identical,  it  should  have  the  same  subjective 
interpretation  as  the  original  signal. 

How  can  we  evaluate  the  quality  of  an  intermediate  representation?  The  usual  method  is  to 
show  that  algorithms  making  use  of  the  representation  arc  effective.  While  it  is  difficult  to  fault 
success,  this  approach  is  limited  by  a  lack  of  generality.  That  is,  the  range  of  algorithms  that  can 
be  demonstrated  is  necessarily  limited,  and  the  reasons  for  their  effectiveness  and  limitations  are 
always  varied. 

Alternatively,  we  can  study  a  representation  mathematically,  to  understand  (1)  the  depen- 
dence of  the  representation  on  the  initial  data,  (2)  the  fibers  of  the  representation,  and  (3)  the  sta- 
bility of  the  transformation.  For  the  first  issue,  the  dep»endence,  it  is  desirable  to  show  some  form 
of  continuity,  so  that  small  changes  in  the  original  data  result  in  small  changes  in  the  representa- 
tion. By  "fibers"  as  studied  in  the  second  issue,  we  mean  the  sets  of  signals  that  map  to  a  single 
representation.  If  the  mapping  is  one-to-one,  then  the  fibers  are  singleton  sets,  and  the  representa- 
tion is  complete.  Finally,  stability  of  the  representation  concerns  continuity  of  the  inverse  map. 
We  would  like  to  show  that  small  changes  in  the  representation,  or  inaccuracies  induced  by 
numerical  implementations,  result  in  small  changes  in  the  fibers  mapping  to  the  respective 
representations. 

Mathematical  analysis  is  complicated  by  the  need  to  define  topologies  on  the  spaces.  For 
example,  the  notion  of  a  small  change  in  the  preimage  sets  (the  fibers)  is  not  a  priori  defined.  The 
metrics  should  correspond  to  an  intuitive  notion  of  similarity  in  the  data,  as  observed  by  human 
subjects.  Further,  analysis  of  a  representation  is  especially  difficult  if  there  are  serious  nonlineari- 
ties  present  in  the  transformation.  The  intent  of  the  results  presented  in  later  sections  of  this  work 
is  to  contribute  to  the  mathematical  analysis  of  representations  involving  zero-crossings  in  scale- 
space. 

As  a  side  benefit  of  a  mathematical  analysis  of  a  representation,  it  is  often  possible  to  devise 
a  reconstruction  (or  approximate  reconstruction)  algorithm.  A  reconstruction  scheme  is  a  right 
inverse,  in  the  sense  that  the  transformation  applied  to  a  reconstruction  should  yield  the  same 
representation  to  which  the  reconstruction  was  applied  (see  Figure  1).  Thus  the  purpose  of  the 
reconstruction  is  to  construct  an  element  in  the  representation's  fiber.  A  reconstruction  can  be 
used  in  a  subjective  or  psychophysical  evaluation  of  a  representation.  Specifically,  a  comparison 
can  be  made  between  an  original  signal  and  a  reconstruction  of  its  representation.  Two  outcomes 
are  possible.  If  one  discovers  an  instance  when  the  reconstruction  and  the  original  are  visually 
distinguishable,  then  the  representation  does  not  encode  the  perceptual  content  of  the  scene,  and 
fails  to  encode  the  differences  that  are  observed.  In  particular,  any  recognition  system  based  on 
the  representation  will  respond  identically  to  the  different  scenes.  Alternatively,  we  may  discover 
that  the  reconstructions  are  generally  indistinguishable  from  the  original,  at  least  in  terms  of  the 
content  for  recognizability  of  the  objects  or  for  the  intended  tasks.  Note  that  it  is  not  necessarily 
expected  that  the  reconstructions  are  numerically  identical.  Indeed,  if  the  reconstructions  are 
numerically  different  but  perceptually  the  same,  then  this  behavior  forms  strong  support  for  the 
utility  of  the  representation  for  building  a  vision  system  that  mimics  some  portion  of  human 
visual  capabilities.  That  is,  patterns  of  indistinguishability  are  evidence  that  the  representation 
carries  information  essential  to  the  interpretation  of  the  associated  signals. 

Ideally,  one  would  like  to  di.scover  and  computationally  simulate  the  same  representation  as 
is  used  by  the  human  system.  Unfortunately,  despite  extensive  study  by  neuroscientists  of  the 
visual  pathway,  complete  computational  models  are  still  elusive.  Suppose  we  were  to  find  a 
representation  that  behaves  well  mathematically,  and  yields  patterns  of  indistinguishability  simi- 
lar to  those  of  the  human  visual  system.  This  would  then  constitute  evidence,  but  not  conclusive 
evidence,  that  the  representation  is  used  by  the  human  visual  system.  Given  the  variety  of  possi- 
ble features  that  could  be  used  in  a  rcprc.scntation,  wc  would  have  to  be  ver\'  lucky  to  pick  pre- 
cisely    the     correct     representation     without     considerable     physiological     evidence.      For 
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Figure  1.  Although  a  good  representation  need  not  be  complete,  it  should  be  a  right  inverse,  in 
that  the  representation  of  a  reconstruction  should  be  the  same  as  the  original  representation. 


representations  based  on  zero-crossings,  to  be  considered  later  in  this  paper,  there  is  some  minor 
evidence  that  there  might  be  cortical  cells  that  respond  to  zero-crossings  [3],  and  considerable 
evidence  (particularly  psychophysical  evidence)  that  zero-crossings  are  not  used  by  the  human 
visual  system.  Accordingly,  our  results,  which  include  instances  of  stable  reconstructions,  pro- 
vide interesting  infomiation  about  potential  intermediate  representations  for  machine  vision  sys- 
tems, but  can  not  be  used  to  make  conclusions  about  the  computations  of  the  human  visual  sys- 
tem. 

1.2.  Scale-space 

Early  work  in  image  processing  and  signal  pattern  recognition  emphasized  the  development 
of  simple  local  operators  for  the  detection  of  features.  For  example,  dozens  of  edge  detectors 
have  been  developed  for  image  analysis  apphcations.  Many  of  the  operators  assumed  a  certain 
resolution  level,  and  apply  the  of)erations  to  a  neighborhood  of  samples  that  is  independent  of  the 
data.  Thus  many  edge  detectors  assume  that  the  edge  can  be  discerned  within  a  three  by  three 
pixel  neighborhood.  Such  operators  frequently  fail  to  find  features  when  the  data  is  sampled  too 
finely,  and  sometimes  fail  when  the  resolution  is  too  coarse.  Nonetheless,  in  controlled  cases, 
feature  operators  can  perform  useful  preprocessing  for  the  analysis  of  signal  data. 

In  order  to  build  into  the  process  a  certain  degree  of  scale  invariance,  and  to  make  the 
feature  detection  operate  dynamically  on  the  data,  the  use  of  multiresolution  data  structures  has 
become  an  important  aspect  of  signal  processing.  One  of  the  earliest  uses  of  multiresolution 
methods  in  computer  vision  was  described  by  Rosenfeld  and  Thurston  [4].  A  simple  construction 
is  to  sample  the  image  at  multiple  resolution  levels,  forming  a  "pyramid"  of  images.  Each  layer 
represents  a  different  scale  for  the  same  image.  Algorithms  are  Uien  developed  that  incorporate 
feature  detector  responses  obtained  by  applying  the  same  operation  to  each  level  of  the  pyramid. 
From  a  single  high  resolution  image,  this  pyramid  can  be  obtained  numerically  by  successively 
blurring  and  subsampling  the  original  image.  Each  operation  of  blurring  and  subsamphng  pro- 
duces an  anti-aliased  version  of  the  same  image  at  a  coarser  resolution.  The  resulting  structure  is 
called  a  "Gaussian  pyramid."  Computational  methods  for  building  Gaussian  pyramids,  and 
other  pyramids,  are  described  by  Burt  and  Adelson  [5]. 

We  will  describe  a  continuous  formulation  of  these  multiresolution  pyramids,  in  order  to 
assist  in  the  mathematical  analysis,  and  to  formulate  scale-space  precisely.  We  begin  with  data 
f(x),  where  xeR".  For  n=\,  fix)  could  be  an  acoustical  signal;  for  n=2,  we  usually  regard /(x) 
as  an  image.  The  case  n=3  arises  with  seismographic  data,  time-varying  imagery,  or  tomography 
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data.  Higher  dimensional  domains  for  the  data  are  possible.  Scale-space  refers  to  a  domain  {x,t) 
for  a  set  of  data,  u  {x,t),  parameterized  by  a  variable  t&  R,  giving  variations  on  the  data  /  {x).  In 
particular,  we  will  always  assume  that  f>0,  and  that  u{x,  0)=f(x).  Further,  uix,t)  will  be  con- 
tinuous in  r,  so  that  for  any  to>0,  uixjo)  gives  a  variant  of /(x).  The  idea  of  scale-space, 
parameterized  by  a  continuous  variable  obtained  as  the  standard  deviation  of  a  Gaussian,  is  due  to 
Witkin  [6]. 

The  natural  framework  for  the  analysis  of  scale-space  formulations  of  multiresolution 
representation  is  in  terms  of  the  heat  equation  [7, 8].  We  define  u  (x.O  to  be  a  bounded  solution  to 
the  heat  equation: 

—  =  Am,  (Heat  Equation) 

u(x,0)=f(x). 

The  solution  is  given  by  convolution  against  the  fundamental  solution  to  the  Heat  Equation, 
which  for  the  domain  R"  is  given  by 

u(x,t)=  j  K{x-y,t)f(y)dy, 

R" 


where 


K(x,t)  =  (4nt)-"'^e-^ ''''*'. 


We  see  that  uix,t)  is  obtained  by  blurring  f(x)  by  increasingly  diffuse  Gaussians,  parameterized 
by  f>0,  with  standard  deviations  o  satisfying  2a^  =4f.  In  computer  vision,  scale-space  some- 
times refers  to  the  (x,a)  variables  that  can  be  used  to  rcparameterizc  the  domain  of  u.  Wc  retain 
the  (x,t)  parameterization  to  keep  the  linear  Heat  Equation  relation  for  the  function  u. 

Convolution  by  Gaussians  is  considered  special  for  many  reasons  [7,9,10].  We  see  from 
the  above  analysis  a  relationship  between  Gaussian  convolution,  the  Heat  Equation,  and  the 
Laplacian  operator.  Of  course,  Gaussian  convolution  enjoys  other  properties;  for  example,  the 
central  limit  theorem  implies  that  Gaussian  convolution  is  easy  to  implement  by  an  iterative  pro- 
cedure. 

In  computer  vision,  in  an  idea  that  dates  back  to  1955  [11],  the  image  data  f(x)  is  often 
filtered  with  the  Laplacian  of  a  Gaussian,  instead  of  filtering  by  a  Gaussian.  Filtering  by  the 
Laplacian  of  a  Gaussian  can  be  written  in  three  ways: 

AK*f=K*Af  =  A{K*f). 

If  we  denote  the  result  by  v(x,t),  we  see  that 

(1)  v(x,t)  is  the  fix)  data  filtered  by  the  Laplacian  of  a  Guassian. 

(2)  v(x,t)  is  the  solution  to  the  Heat  Equation  with  initial  data  A/. 

(3)  vix,t)  is  Au(jc,f),  where  u  is  the  solution  to  the  Heat  Equation  with  initial  data/(x). 

We  will  refer  to  v  as  the  scale-space  function  of  the  data  /. 

A  difference  of  two  Gaussians  is  often  used  in  place  of  the  Laplacian-of-Gaussian.  We  see 
from  the  Heat  Equation  formulation  the  basis  for  the  approximation.  Namely,  since  K{x,t)  is 
itself  a  solution  to  the  Heal  Equation, 

(AAf)  (x,t)  =  ^{x,t)  =  \\m(Kix,t+z)  -  K{,x,t))lx. 

That  is,  the  difference  of  Gaussians  is  a  good  approximation  to  AA"  as  the  separation  between  the 
spread  of  the  two  Gau.ssians  approaches  zero  (and  tlic  difference  is  scaled).   In  actual  u.se,  it  is 
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more  frequent  to  use  a  difference  of  Gaussians  where  the  ratio  of  the  two  scales  in  the  a-variable 
is  given  by  1.6,  but  the  rationale  for  this  choice  is  not  the  degree  of  approximation  to  the 
Laplacian-of-Gaussian  [12]. 

1.3.  Zero-crossings 

The  zero  set  of  v(x,t)  is  the  point  set  in  (x,t)  where  v  =  0.  The  set  might  be  empty  (for 
instance,  if  /  is  subharmonic  or  superharmonic;  [13])  everything  (if /is  harmonic;  [14]),  or  a 
proper  subset  of  (x,t)  space.  In  the  latter  case,  zeros  can  be  isolated  points,  lines,  and  surfaces 
(but  never  regions).  We  distinguish  components  of  the  zero  set  which  form  manifolds  of  codi- 
mension  one: 

Definition:  The  zero-crossings  ofv(x,t)  refers  to  the  point  set 

d{(x,t)\v(x,t)<0}  n  d{(x,t)\vix,t)>0].  U 

Zero-crossings  have  been  suggested  for  segmentation  of  imagery  by  edge  detection  [12],  and  for 
stereo  matching  and  motion  correspondence  between  pairs  of  images  (e.g.,  [15]).  It  has  also  been 
suggested  [10]  that  the  zero-crossings  are  a  nearly  complete  representation  of  A/.  Marr  also  sug- 
gested that  the  zero-crossings  might  form  a  complete  representation  [1],  but  also  thought  that 
perhaps  the  zero-crossing  information  would  need  to  be  supplemented  with  gradient  data  along 
the  zero-crossings.  These  suggestions  form  the  basis  of  the  interest  in  zero-crossings  as  an  inter- 
mediate representation  for  images. 

It  should  be  noted  that  there  is  a  long  and  distinguished  history  of  the  analysis  of  zero- 
crossings  in  mathematics.  Kedem  supplies  a  survey  together  with  some  advances  [16].  In  this 
and  related  work,  analysis  is  generally  restricted  to  signals  of  one  variable,  and  the  representation 
is  based  on  the  local  statistics  of  zero-crossings  (the  number  of  crossings  per  unit  length).  Issues 
such  as  the  local  density  of  zero-crossings  become  very  important  when  zero-crossings  are  used 
for  stereopsis  computation  with  image  pairs.  Presumably,  this  body  of  mathematical  work  is 
quite  relevant  to  image  representation,  but  has  not  been  applied  to  computer  vision  in  a  major 
way. 

Witkin  [6]  observes  that  zero-crossings  in  scale-space  evolve  as  t  increases,  and  are  never 
created  at  some  nonzero  t .  This  property,  discussed  in  [9]  and  in  [10],  ensures  that  zero-crossing 
surfaces  are  nested,  one  within  another,  enclosing  regions  containing  the  face  (r  =  0],  or  forming 
a  sheet  meeting  the  face  {t  =0}  and  extending  to  t  =  °°.  The  property  can  be  given  a  precise 
statement: 

Evolution  property  of  zero-crossings:  Let  C  be  a  connected  component  of  the  set  of  zero- 
crossings  in  the  domain  [(x,t)\xG^'',  Ti<t<T2},  where  0<r]<r2.  Then 
C  n{(x,t)\t  =  Ti]:A0.  U 

The  practical  import  of  this  property  is  that  the  zero-crossing  information  can  be  simplified, 
and  represented  in  a  symbolic,  or  approximate  way,  by  describing  the  nesting  and  some  other 
simple  features  of  the  zero-crossing  surfaces.  The  hopt  is  that  the  simplified,  or  symbolic 
representation  of  the  zero-crossing  surfaces  will  also  suffice  to  form  a  complete  representation  of 
the  original  data,  at  least  for  image  analysis  purposes.  For  instance,  Johansen  considers  the 
representation  of  signals  by  the  "toppoints"  of  the  zero-crossings  in  scale-space  (with  t  extended 
to  -co)  [17].  To  date,  the  hope  of  simplifying  the  representation  of  the  fingerprint  of  zero- 
crossings  has  borne  little  practical  fruit,  but  has  nonetheless  sparked  a  wide  range  of  interesting 
mathematical  analysis. 

In  the  next  section,  we  consider  some  of  the  mathematical  results  related  to  the  evolution 
property  and  the  issue  of  completeness  of  the  representation  by  zero-crossings. 
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2.  Analytical  results 

2.1.  Evolution  property 

In  this  section,  wc  show  that  the  evolution  property,  applied  to  level-crossings  (and  not  just 
zero-crossings)  is  equivalent  to  the  classical  Maximum  Principle  for  parabolic  partial  differential 
equations.  That  is,  we  will  show  that  the  evolution  property  can  be  derived  simply  from  the  max- 
imum principle,  and  moreover,  that  any  scale-space  construction  that  obeys  the  evolution  pro- 
perty for  all  level-crossings,  subject  to  minor  restrictions,  will  satisfy  a  maximum  principle. 

The  classical  maximum  principle  for  the  solution  to  the  parabolic  Heat  equation  du/dt  =  Au 
states  (see,  e.g.,  [18-20] ): 

Maximum  Principle:  Let  DqR"  be  open  and  bounded.  Suppose  «  is  a  solution  to  the  Heat 
equation  in  7=  {(x,t)\xeD,0<t<T}  of  class  C^  and  is  continuous  in  the  closure  T.  Then  u 
assumes  its  maximum  at  some  point  (x,t)  for  which  either  j:e  3D  or  /  =  0. 

The  maximum  principle  holds  for  more  general  parabolic  partial  differential  operators,  but 
certainly  not  for  all  parabolic  equations.  It  is  not  known  (to  the  authors'  knowledge)  how  to 
characterize  parabolic  equations  giving  a  maximum  principle.  For  example,  the  maximum  prin- 
ciple does  not  hold  when  the  biharmonic  operator  is  used:  du/dt  =  A^u.  Generally,  the  elliptic 
operator  will  be  second  order  in  order  for  a  maximum  principle  to  hold. 

We  will  assume  that  whenever  a  maximum  principle  holds  a  similar  minimum  principle 
also  is  given.  For  linear  parabolic  operators,  such  as  the  Heat  equation,  a  minimum  principle  fol- 
lows from  the  maximum  principle.  However,  since  we  can  envision  more  general  operators,  we 
will  henceforth  assume  that  the  statement  that  a  maximum  principle  holds  means  that  both  a  max- 
imum and  minimum  principle  hold. 

Next,  let  us  denote  the  construction  of  the  representation  of  the  image  /  in  scale-space  by 
the  operator  v  =  5/  In  the  previous  section,  we  described  the  construction  of  v,  given  /e  L"'(R'), 
according  to 

v(x,t)=  j  AK(y-x,t)f(y)dy  . 

However,  we  now  wish  to  permit  more  general  constructions,  allowing,  for  example,  /  to  be 
defined  on  an  irregular  or  bounded  domain,  or  allowing  constructions  involving  differential 
operators  with  nonconstant  coefficients.  However,  we  will  need  to  place  restrictions  on  the 
growth  of  functions  v=5/as  Ix  I— »«  and  on  the  boundary  of  the  domains.  Specifically,  we  will 
assume  that  the  domain  of /is  D  (which  may  or  may  not  be  bounded),  and  that  the  corresponding 
scale-space  will  have  (x,t)eDx[0,°o),  There  is  a  class  of  permissible  images  /for  which  the 
scale-space  operator  v  =  S/is  defined  and  we  assume  that  for  v's  constructed  in  this  fashion; 

(/)    vgC(Dx(0,«)) 

(a)    v(x,f)  =  0for;ce9D,  f>0; 

(in)    v(x,t)  ^  0  as  I  j:  I  — »  0°  uniformly  for  r  >  0  . 

Note  that  this  means  that  for  the  Laplacian-of-Gaussian  construction  given  above,  we  must  limit 
application  to  functions /that  tend  to  zero  as  U  I  — »  «>.  The  space /eL^(R^)r\L~(R*)  suffices, 
for  example.  This  restriction  is  stronger  than  necessary,  but  not  terribly  objectionable. 

Our  main  result  in  this  section  is  that  a  scale-space  construction  v=S/.  given  the  boundary 
restrictions,  satisfies  the  evolution  property  for  level-crossings  if  and  only  if  a  maximum  principle 
holds  (meaning  both  a  maximum  and  minimum  principle). 

The  result  is: 
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Theorem  1:  Consider  a  scale-space  construction  v=S/ satisfying  (/),  (ii)  and  070  above  for  all 
admissible  functions  /.  Then  the  following  are  equivalent: 

(1)  The  maximum  principle  holds  for  functions  v  =5/ obtained  from  admissible  /s. 

(2)  The  evolution  property  holds  for  level-crossings  of  solutions  v  =S/ obtained  from  admis- 
sible/s. 

Proof:  We  first  show  that  the  maximum  principle  implies  the  evolution  property.  For  if  the  evo- 
lution property  fails  for  a  level-crossing  /  of  vix,t),  v  =  Sf,  then  there  is  a  component  C,  and 
bounds  ti,t2,  where  C  is  a  connected  component  of  {v(x,t)>  l\t]<t<t2)  or 
[v(x,t)  <  1 1  f  i<r<r2)  such  that  C  does  not  meet  the  plane  {f=f  i ).  If  1*0,  then  the  component  C 
must  be  bounded,  since  Iv  1-^  uniformly  as  \x  1^°°.  Thus  the  supremum  (or  infemum)  of 
v(x,t)  is  attained  at  a  finite  point  (xo.ro)  in  C.  If  /=0,  then  the  supremum  (or  infemum)  in  C  is 
nonzero,  and  once  again  attained  at  a  point  (xo.^o)  of  C.  Since  (xo,fo)  cannot  lie  on  the  level- 
crossing  surface  (assuming  v  is  nonconstant),  it  must  either  be  in  the  interior  of  C,  or  lie  on  the 
top  surface  [1=12}  within  C.  In  either  case,  we  can  place  a  bounded  cylinder  within  C  with 
(xo,to)  in  the  interior  or  top,  in  violation  of  the  maximum  principle. 

To  show  that  the  evolution  property  for  level-crossings  yields  the  maximum  principle,  an 
even  simpler  argument  suffices.  For,  suppose  that  Dx[ti,t2]  is  a  bounded  cylinder,  with  a  max- 
imum (or  minimum)  at  a  point  (;co,ro)  in  the  interior  or  top  of  the  cylinder,  in  violation  of  the 
maximum  principle.  Either  way,  there  is  a  value  /  less  than  the  maximum  (or  greater  than  the 
minimum)  at  (xo.To),  but  greater  (or  less  than)  the  values  of  v  on  the  sides  and  bottom  of  the 
cylinder.  Thus  a  component  of  the  level-/  crossing  within  the  wedge  t\<t<t2  lies  entirely  within 
the  interior  of  the  cylinder,  and  thus  does  not  meet  the  plane  {t=ti).  Thus  the  evolution  property 
fails  if  the  maximum  principle  is  violated.  ■ 

As  a  result  of  the  theorem,  we  have  an  evolution  property  for  scale-space  constructions 
more  general  than  Laplacian-of-Gaussian  filtering.  For  example,  let  D  be  a  bounded  domain,  and 
consider  scale-space  T  =  Dx[0,«')  with  functions  constructed  as  follows: 

For /g  L"(D) ,  solve  for  mg  C(T), 

Am  =  4^    inDx(0,oo) 
at 

uix,0)=f(x). 
u(,x,t)=f(x),xedD,  t>0. 
Then  let 

v(x,t)=^(x,t). 

We  observe  that  v=0  on  the  sides  of  the  scale-space  cylinder,  and  that  the  maximum  principle 
holds  since  v  satisfies  the  Heat  equation.  Thus  level-crossings  of  v  evolve  as  r  increases,  and 
always  meet  the  base  {f=0). 

Under  more  restrictive  assumptions,  one  can  show  that  the  evolution  property  requires 
Laplacian-of-Gaussian  scale-space  construction  [9].  The  assumptions,  however,  require  that  the 
domain  be  all  of  Euclidean  space,  and  that  the  underiying  equation  have  constant  coefficients. 
Further,  proofs  of  the  evolution  property  have  frequently  assumed  that  a  level-crossing  com- 
ponent not  meeting  {f=0}  will  have  an  extremal  lower  point,  which  is  an  invalid  restriction. 
Given  the  equivalence  of  the  evolution  property  to  the  maximum  principle,  any  proof  that  does 
not  cite  the  maximum  principle  or  essentially  redo  the  proof  is  suspect.  Since  the  proof  of  the 
maximum  principle  is  slightly  delicate,  especially  in  the  absence  of  strong  regularity  assump- 
tions, the  former  course  seems  more  appropriate. 
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A  version  of  the  evolution  prop)erty  is  easy  to  establish  when  the  scale-space  is  a  discrete 
domain.  For  simplicity,  let  us  treat  the  case  of  one  space  dimension,  and  assume  that  we  are 
given  data  fi,  i  =  •  •  •  ,-1,0,1,  ■  ■  • .  We  obtain  scale-space  data  v,jt,  where  k  is  the  scale,  or 
level,  k  =  0,],  ■  ■  ■  .  Tlie  method  by  which  v  is  constructed  is  unimportant  at  this  point;  a  max- 
imum principle  states  that  in  any  rectangular  grid  of  lattice  pwints  [(i,k)\ii<i<i2,ki<k<k2],  the 
maximum  of  v,  ^  must  occur  on  either  the  bottom  k=ki  or  the  sides,  i=ii  or  i=i2-  The  maximum 
principle  will  surely  hold  true  if  v  satisfies  a  construction 

Vi,*+i  =a,-i,tv,_i_i  -I-  a,.*v,_i  +  a,+i_tv,+]jt . 

where 

a,-\.k  +  a,_k  +  a,+i,i  =  1     and  all  a,,*  >  0. 

Given  the  maximum  principle  for  discrete  data,  we  can  easily  obtain  an  evolution  property.  The 
evolution  property  will  state  that  for  any  4-connecied  component  of  grid  points  of 
{0.^)1  V,  *  >0,  *i  <k<k2  ),  this  component  will  include  points  on  the  bottom  level  k=ki. 

2.2.  Completeness  when  gradients  are  included 

We  now  show  that  the  zero-crossings,  when  supplemented  with  gradient  data  along  the 
zero-crossing  boundaries,  are  sufficient  to  reconstruct  the  original  function  fix).  Actually,  the 
method  is  used  to  reconstruct  Af(x).  However,  if  we  assume  that  / (x)-*0  as  I x  I  — >«,  then  / (x) 
can  be  reconstructed  by  solving  Poisson's  equation.  Details  of  these  ideas  were  reported  earlier 
in  a  technical  report  [21].  In  the  same  report,  it  is  observed  that  completeness  of  the  representa- 
tion involving  gradients  along  the  zero-crossings  is  easily  established,  by  a  non-constructive 
proof,  using  the  Hopf  maximum  principle  for  the  Heat  Equation.  We  instead  provide  a  construc- 
tive proof  below,  although  it  should  be  noted  that  the  constructive  method  is  unstable. 

The  use  of  gradient  data  for  the  representation  also  appears  in  [7],  but  the  gradient  data 
there  is  not  limited  to  the  zero-crossings.  The  use  of  gradient  data  along  zero-crossings  is  dis- 
cussed in  [22].  Many  researchers  have  noted  from  a  casual  observation  of  zero-crossings  of 
image  data  that  zero-crossings  with  large  gradient  magnitudes  are  of  greater  significance  than 
those  with  low  gradient  magnitudes. 

2.2.1.  Continuous  Case 

Specifically,  let  Q  be  a  bounded  connected  component  of  [(x,t)  I  f>0,  v{x,t)^},  and 
denote  by  D  the  set  {x&  R"  I  (x,  0)e  Q],  and  by  F  the  zero-crossing  3Qn{r>0}.  Let  t  be  a  value 
such  that  x>sup(f  l(.r,r)€Q).  Next,  denoting  Af(x)  by  g(x),  we  set  g(x)  =  g(x)  for  xgD,  and 
^(a:)  =  0  elsewhere.  Let  ^U)  be  the  |(.x)  data  blurred  to  the  level  t: 

bix)=  JKix-y,T:)g(y)dy. 

R" 

Finally,  for  fixed  y  and  x  define 

w(x,t)  =  K(x-y.z-0. 

Theorem  2.  Suppose  that  the  data  g  (x),  xe  D,  is  contained  under  a  bounded  zero-crossing  sur- 
face r  of  the  scale-space  function  v.  Then  b(x),  and  hence  g{x)  for  xeD,  can  be  reconstructed 
from  the  data  F  and  Vv  along  F. 

Proof.  Observe  that 

dt 
Using  this  fact  and  Green's  theorem  we  have 
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3^        A 
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■vdxdt 


■w  dxdt  +  j 
=  I  [(vH')^  +  (vVH'-wVv)-n]do 

where  (n,0(j:,r)  is  the  outward  pointing  surface  normal  to  r  at  (x,r)  and  da  is  surface  area  meas- 
ure. The  gradient  V  is  with  respect  to  spatial  coordinates  only.  Using  the  facts  that 
(n,Q  =  (0,-1)  at  r  =  0  and  v  =  0  on  F,  the  equation  reduces  to 

Jv(;c,  Q)w{x,  G)dx=-\w{x,t)(yvn)ix,t)da. 
D  r 

Using  the  definitions  of  |,  b  and  w  above,  the  equation  becomes 

b{y)=  \K{x-y,x-t){yvn){x,t)do. 
r 

Thus  given  the  zero-crossing  T  and  ^v{x,t)  for  (jc,f)e  T,  the  blurred  data  b(x)  can  be  con- 
structed by  a  simple  linear  process.  The  original  data  g{x)  can  be  reconstructed  by  deblurring  the 
ft(;c)data[23].  ■ 

Deblurring  is,  of  course,  a  classic  unstable  process.  The  situation  is  not  hopeless,  however, 
since  g{x)  has  known  compact  support,  which  might  be  used  to  advantage,  and  also  since  errors 
that  occur  are  predominantly  in  high  frequency  components,  which  might  not  be  as  essential  to 
visual  interpretability.  Nonetheless,  the  method  is  still  unstable. 

However,  the  foregoing  analysis  applies  only  to  the  case  where  fl  is  bounded.  Next  con- 
sider an  unbounded  component  i^  of  {x,t)  space  with  the  zero-crossings  removed.  Let  x  be 
greater  than  the  range  of  t  in  aU  of  the  bounded  components  of  that  space.  Finally,  let 
Q^=Q.n  [{x,t)\  0<r<T}.  We  repeat  the  previous  argument,  replacing  integration  over  Q.  by 
integration  over  Q.^.  There  is  now  in  addition  to  D  =  [xe  R"  t  {x,  0)€  Qx  1 .  another  boundary  por- 
tion D  X  =  {A:e  R"  I  (;c,  t)e  fix }  to  consider  in  the  surface  integral.  Thus  we  obtain 

\v{x,Qi)w{x,Q)dx+     \     w{x,t)(yvn){x,t)d<5=y^''^^   ^  ^ 

D  rn(«T)  [0  otherwise. 

Here  we  have  used  the  fact  that  w(y,x)=K (x-y,0)  is  a  delta  function  centered  at  x=y.  With  g(x) 
and  bix)  defined  as  before,  we  have 

b(y)=      J     K(,x-y,x-t)iVvn)(x,t)da,  (y,x)4D^. 

rn(/<T) 

Since  v(y,x)  is  imknown,  the  equation  gives  no  useful  information  when  (y,x)eD^.  The  missing 
data  could  be  supplied  by  analytic  continuation  of  ^Cy)  from  (y,x)4D^,  assisted  by  the  fact  that 
the  above  equation  holds  with  the  unknown  but  (hopefully)  small  error  v(y,x)  when  (y,T)€Dx. 
Once  b(y)  has  been  extended  in  this  way,  g(x)  is  once  again  found  by  deblurring. 

The  lesson  of  this  section,  ultimately,  is  that  even  for  bounded  zero-crossings  supplemented 
with  gradient  data  along  the  zero-crossing,  reconstruction  by  the  indicated  method  is  still 
unstable. 

2.2.2.  Discrete  Data 

It  is  interesting  that  the  analytical  result  above  can  be  converted  to  a  discrete  form,  without 
involving  any  discrete  approximations.  This  is  because  Green's  theorem  can  be  converted  to  a 
discrete  simi.  For  a  signal  or  image  defined  on  a  discrete  lattice,  we  first  define  a  notion  of  a 
Gaussian  pyramid  (actually,  a  "monolith"),  and  then  a  Laplacian  pyramid,  and  then  the  concept 
of  zero-crossings  in  the  discrete  scale-space.   Knowledge  of  the  zero-crossings  together  with 
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gradient  data,  in  the  formulation  below,  will  translate  to  knowledge  of  the  values  of  the  scale- 
space  function  on  a  discrete  set  of  points,  with  the  points  located  on  either  side  of  a  zero-crossing. 
Since  values  are  given  on  both  sides  of  a  zero-crossing,  the  gradient  values  are  implicitly  encoded 
in  the  representation. 

For  simplicity,  we  treat  the  case  of  one  unbounded  space  dimension,  although  the  results 
extend  easily.  We  are  given  data/,,  /=•••  ,-1,0,1,  ••■ ,  and  define 

1  ^  1  ^      1  . 

gi  -  ~7ji-\  -  ~zJi  +  T/i+1  ■ 
4  2  4 

We  define  the  filtered  data  v^^|^  recursively: 

v..o  =  «.. 


1  1 

M  - 

We  also  define  the  blurring  kernel 


4  2  4 


'^'•*  -  ^ 


2k 
J+k. 


,    -k<i<k,   k>0, 


and  AT,  *  =  0  elsewhere,  for  k>0.  Both  v  and  K  satisfy  a  discrete  version  of  the  Heat  Equation, 
namely 

1  1  1 

v,.*+i-v,,i  =  -rv,-i,*  -  — v,.t  +  — v.+i,i. 
4  2  4 

It  is  not  hard  to  prove  a  discrete  analogue  of  the  evolution  property  for  zero-crossings.  The 
key,  as  indicated  in  Section  2.1,  is  a  discrete  version  of  the  maximum  principle,  which  is  easy  to 
establish. 

To  formulate  the  reconstruction  method,  let  ti  be  a  bounded  4-connectcd  collection  of  pix- 
els (/,/(:)  with  a  nonempty  set  D  ={;  I  (j,  0)eQ).  Let  r  be  an  upper  bound  r> max {/k  I  (/,i)en), 
and  define  i>,  to  be  the  data  gi,ieD,  blurred  to  level  T: 

Finally,  denote  the  set  of  pixels  on  the  boundary  of  Q  by 

3(±i,o)"={('.*)en  I  (/±l,^)^n 

3(0.1)"  ={(«.'t)en  I  O'./t-t-D^n}. 

5(o.-i)i^  ={('.*)€ i2  I  *>0.  (j,A:-iyn}. 
Then  simple  but  messy  algebra  allows  us  to  show 
Theorem  3. 


4^=  I     I 

£=-I.l     (l,*)€8(£.0)n 


V,,i+V,>5.,t 


(Ki+t-j.T-k-\  - ^i -j.T-k-l ) 
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f^i+t-jJ-k-l+^i-j.T-k-\  ,  . 
T (v,+e.*  -  V,-,*) 

+     Z     '*^<.*+i  f^i-jj-k-i 

(i.*)€  3(0,1)0 

(i.*)e  3(0,-1)0 

To  reconstruct  data  by  the  above  equation,  choose  a  connected  component  of 
{(/,/t)  I  v(i,k)>0],  (or  respectively  <0  ).  If  the  component  extends  to  infinity  in  either  coordi- 
nate, truncate  the  domain  to  become  a  convenient  bounded  collection  of  pixels,  and  denote  the 
result  by  Q.  We  store  the  sets  3(±i,o)^.  3(o.±i)^.  and  D  as  defined  above.  For  pixels  (i,k)  in 
9(i_o)J^  (respectively  3(-i,o)f^  ).  we  store  the  information  v,  *  and  v,+,jk,  (respectively  v,,*  and 
v,-i,*).  For  pixels  0',^:)  in  3(o.i)^,  we  store  the  data  v,  ^+1  and  for  3(o,-i)^  pixels  we  store  v,jt. 
Using  the  above  equation,  we  choose  a  T  and  reconstruct  the  blurred  data  bj.  To  reconstruct  the 
data  gi  for  j'eD,  it  suffices  to  deblur  the  b,  data  by  solving  for  gi  in  the  linear  equations  defining 
bi.  In  fact,  the  system  is  overdetermined,  although  still  poorly  conditioned,  especially  if  ID  I  or 
T  is  large. 

In  order  to  make  the  computations  feasible,  it  is  necessary  to  modify  the  formulas  for  a 
bounded  spatial  domain.  For  example,  a  common  approach  to  constructing  a  discrete  scale-space 
function  involves  solving  a  bounded  domain  problem,  with  -N<i<N,  by  setting  v,  *  =  0  for 
i  =  ±A^.  The  blurring  kernel  K  is  changed  by  this  modification,  but  the  proposition  carries  over 
with  little  change. 

There  is  another  version  of  the  formula  for  bj  above,  and  thus  for  the  reconstruction  of  the 
image  data,  where  the  stored  information  consists  of  the  precise  locations  of  the  zero-crossings, 
assuming  a  linear  interpolation  of  values  between  grid  points.  The  gradient  data,  in  the  fonm  of 
the  difference  of  the  neighboring  grid  values,  still  needs  to  be  included.  With  this  infonmation, 
the  fomiulas  simplify  some,  but  the  information  in  the  representation  is  clearly  equivalent,  so  we 
omit  the  details. 

3.  Completeness,  and  a  brief  survey  of  reconstructions 

Under  certain  restrictions  on  the  class  of  signals,  it  can  be  shown  that  the  zero-crossings  in 
scale-space  form  a  complete  representation  of  the  signal  data.  There  have  been  many  different 
studies  and  theorems  along  these  lines  for  applications  to  computer  vision.  An  oft-cited  study, 
for  the  case  of  signals  in  one  dimension,  is  provided  by  Logan  [24].  Logan's  theorem  requires 
that  the  function  be  bandpass,  and  satisfy  certain  technical  conditions  concerning  the  signal's  Hil- 
bert  transform.  Further,  the  result  is  intended  for  the  reconstruction  of  a  function  represented  by 
the  zero-crossings  at  exactly  one  level  of  resolution.  While  feasible  as  a  complex  analytic  result, 
our  interest  here  is  in  the  more  practical  possibility  of  reconstruction  given  zero-crossings  at  mul- 
tiple scales  of  resolution  in  scale-space. 

More  complex  analytic  results  can  be  used  when  the  image  data  is  polynomial.  In  fact, 
when/(;c)  and  thus  Af(x)  is  a  polynomial  in  xsH",  then  v(x,t)  is  a  polynomial  in  (x,t)e¥L"*\ 
Accordingly,  the  zero-crossings  are  part  of  the  analytic  varieties  of  the  polynomial  v  as  studied  in 
algebraic  geometry.  It  is  well  known  that  the  varieties  in  C  determine  the  complex  polynomial 
defined  on  n  complex  variables.  It  is  not  as  commonly  used,  but- nonetheless  true,  that  an  n- 
dimensional  subportion  of  the  intersection  of  the  analytic  variety  with  R"""'  detenmines  an  irredu- 
cible polynomial.  A  proof  of  this  result,  supplied  by  Professor  D.  Mumford,  is  given  in  an  appen- 
dix. Related  proofs  are  given  by  Huang  and  Sanz  [25]  and  Curtis,  Oppenheim,  and  Lim  [26]. 
Thus  the  case  of  polynomial  data  can  be  setUed  with  algebraic  geometry. 
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However,  these  results  depend  heavily  on  the  assumption  that  the  data  is  polynomial.  It  is 
obvious  that  there  are  many  different  functions,  even  smooth  functions,  that  have  identical  zero- 
crossings.  Interestingly,  the  completeness  in  the  case  of  polynomial  data  applies  equally  well  to 
the  scale-space  function  v{x,t)  as  to  the  original  data/(j:),  or,  for  that  matter,  any  single  level  of 
the  scale-space  function  vixjo)-  That 's.  the  polynomial  /(j:)  is  determined  by  its  zero-crossings 
(providing  it  is  irreducible  as  a  polynomial),  a  level  of  the  scale-space  function  v(j:,fo)  is  deter- 
mined by  its  zero-crossings,  providing  it  is  irreducible,  or  the  scale-space  function  v(x,f)  is  deter- 
mined by  its  zero-crossings,  providing  it  is  irreducible.  In  any  of  the  above  cases,  f(x)  can 
theoretically  be  reconstructed.  However,  it  is  in  some  sense  more  likely  that  v(x,t)  will  be 
irreducible,  and  we  expect  the  reconstruction  of  /  from  the  function  v  to  be  more  stable  than 
reconstruction  from  a  single  level  of  v,  or  from  the  zero-crossings  of/ alone. 

Yuille  and  Poggio  [27]  prove  a  similar  result:  showing  that  when  Af(x)  and  hence  f(x)  is 
polynomial  in  x,  and  if  n=\,  then  reconstruction  from  zero-crossings  is  theoretically  possible. 
They  refer  to  the  validity  of  the  observation  for  larger  n.  Their  method  relics  on  an  expansion  of 
the  analytic  struction  of  a  zero-crossing  contour,  using  arbitrarily  high  derivatives  of  the  contour, 
at  two  points  of  the  zero-crossing  in  scale-space.  While  the  method  is  constructive,  it  is  not  likely 
to  lead  to  stable  reconstructions. 

Curtis,  Shitz,  and  Oppenheim  [28]  and  Sanz  and  Huang  [25]  extend  the  above  results  to  the 
case  when/(x)  is  a  bandlimited  function  which  is  irreducible  as  an  entire  function.  These  results 
arc  likewise  of  an  algebraic  geometry  nature,  and  rely  on  the  fact  that  there  is  a  finite  Fourier 
series  of  the  bandlimited  function  that  provides  a  polynomial  of  several  complex  variables  for 
representing  /. 

Since  the  determination  of  a  polynomial  by  its  varieties  is  essentially  an  anal>tic  continua- 
tion result,  stability  of  the  reconstruction  is  unlikely.  That  is,  small  errors  in  measurement  of  the 
zero-crossings  could  lead  to  arbitrarily  large  errors  in  the  determination  of  Af(x).  Put  differently, 
there  can  be  widely  different  initial  data  leading  to  neariy  identical  zero-crossing  data. 

Worse,  settling  the  case  for  pwlynomial  data,  or  even  irreducible  entire  bandlimited  data, 
says  little  about  the  general  case  of  continuous  initial  data.  Although  the  Stone-Weierstrass 
theorem  says  that  a  continuous  function  can  be  uniformly  approximated  by  a  polynomial  on  a 
compact  set,  the  zero-crossings  depend  on  the  initial  data  globally,  and  the  dependence  can't  be 
localized.  Further,  the  lack  of  stability  means  that  the  approximation  is  irrelevant.  The  situation 
is  similar  to  the  fact  that  a  polynomial  of  a  single  variable  with  all  real  roots  is  determined  by  its 
zeros,  but  that  given  all  the  zeros  of  a  continuous  function,  one  knows  nothing  more  than  the 
zeros. 

Based  on  these  results  of  completeness  under  certain  restrictions,  there  have  been  a  number 
of  attempts  at  reconstructions  from  zero-crossings.  The  hope  is  that  natural  image  data  will 
sufficiently  restrict  the  class  of  admissible  functions  that  some  stable  reconstruction  method  can 
be  found.  We  briefly  survey  below  three  such  attempts.  In  section  4,  we  provide  a  formulation 
for  a  new  attempt  at  reconstruction. 

3.1.  Curtis-Oppenheim 

Susan  Curtis  and  Alan  Oppenheim,  using  results  that  show  that  bandlimited  images  are 
completely  represented  by  their  zero-crossings  under  certain  conditions,  present  some  experi- 
ments with  reconstruction  [29].  In  their  formulation,  they  begin  with  a  bandlimited  image 
f  {x,y),  and  record  a  large  number  of  locations  in  the  U.v)  domain  where  /crosses  some  thres- 
hold (such  as  zero,  assuming  that  /has  both  positive  and  negative  regions).  These  locations  are 
recorded  with  great  precision,  to  ultra-subpixcl  accuracy.  The  image  is  reconstructed  from  this 
information.  Thus  in  their  work,  the  image  is  recoasiructcd  from  a  number  of  level-crossing  con- 
tours. Clearly,  the  assumption  that  the  image  is  bandlimited  is  crucial,  for  otherwise  there  are  an 


Page  v. 


Hummel  and  Moniot 


infinity  of  images  having  precisely  the  same  level-crossings.  In  essence,  the  work  constitutes 
reconstruction  from  thresholded  imagery,  but  with  the  understanding  that  the  image  is  bandlim- 
ited,  and  that  the  level-crossing  contours  are  recorded  with  great  accuracy. 

The  method  begins  by  writing  the  discrete  image  /  as  the  inverse  Fourier  transform  of  its 
discrete  Fourier  transform  F: 


J    N-\N-\ 

N    „=ov=o 


This  function  is  interpolated  to  an  entire  bandlimited  periodic  function,  /  (,x,y),  given  by 

N-\N-l 


1      /V-lrtf-l 


A'    „=ov=o 

It  is  the  locations  of  zeros  of  f(x,y)  that  are  recorded,  say  (x„,}'„),  n=l,  ■  ■  ,K.  We  can  then 
write  the  K  equations 

S' x'F(a,v)e'"'""'^e'""^""'  =  0,  n  =  1,2,  •  •  •  ,A:. 

u=Ov=0 

If  the  mean  value  (J.  of/  {x,y)  is  not  zero,  then  the  additional  equation  F(0,0)  =  liN^  is  included. 
Otherwise,  one  more  equation  is  needed,  giving  the  location  and  value  of  f(x,y)  at  a  single  loca- 
tion where  /is  nonzero.  Assuming  the  former  case,  we  have  K  linear  equations  in  the  unknowns 
F{u,v),  (m,vMO,0),  where  (m,v)  ranges  over  the  known  (bandlimited)  support  of  F.  By  using 
many  more  equations  than  unknowns,  and  computing  the  least  squares  solution,  reasonable 
reconstructions  are  obtained.  However,  in  the  examples  shown,  the  number  of  independent 
nonzero  spectral  components  is  roughly  200.  Further,  the  locations  of  these  spectral  components 
were  chosen  carefully,  in  advance,  and  the  image  to  be  reconstructed  was  obtained  from  a  full- 
resolution  image  by  setting  all  except  the  200  Fourier  coefficients  to  zero.  By  tiling  Fourier  space 
(as  Rotem  and  Zeevi  do;  see  below),  a  general  image  could  be  reconstructed  by  this  method  using 
the  zero-crossings  of  multiple  images,  each  obtained  by  filtering  the  initial  image.  However, 
since  a  256  by  256  image  contains  64K  independent  spectral  components,  it  is  not  clear  how 
many  levels  will  be  needed.  Nonetheless,  the  method  provides  a  theoretically  possible  means  for 
representation  and  reconstruction  from  multiple  images  of  zero-crossings. 

3.2.  Rotem-Zeevi 

Doron  Rotem  and  Yehoshua  Zeevi,  in  a  paper  in  which  they  extend  the  Logan  theorem  for 
application  to  two-dimensional  images,  experiment  with  reconstructions  from  zero-crossing  data 
[30].  Their  method  is  based  on  first  decomposing  the  given  image  f{x,y)  into  a  sum  of  images 
/*(^.>').  with  each  /t  obtained  from  /by  a  filtering  operation  in  Fourier  space.  Each  /^  has  the 
property  that  its  discrete  Fourier  transform  has  nonzero  coefficients  in  a  region  that  is  bandpass  of 
one  octave  or  less  in  at  least  one  of  the  two  spectral  dimensions.  Further,  the  low-frequency  com- 
ponents are  omitted  from  the  reconstruction  process,  and  simply  stored  as  part  of  the  representa- 
tion, and  an  attempt  is  made  to  reconstruct  only  a  low-pass  version  of/,  where  aU  frequencies 
above  half  of  the  maximum  are  omitted.  Thus  for  a  256  by  256  image,  where  spectral  frequen- 
cies range  from  -128  to  127,  only  the  frequencies  in  the  range  -64  to  64  are  included  in  the 
reconstruction,  and  the  spectral  frequencies  where  both  components  u  and  v  satisfy  ImI<8, 
I V  I  <8,  are  stored  separately,  and  excluded  from  the  reconstruction.  The  remaining  spectral 
domain  is  decomposed  into  19  regions,  so  that  there  result  19  images,  /*,  k=\,  ••  •  ,19.  If /o 
denotes  the  low-pass  image  composed  of  frequencies  below  8,  and  if  g{x,y)  is  the  high-pass 
image  composed  of  frequencies  above  64  (in  either  dimension),  then 
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By  keeping  the  filters  yielding  the  decomposition  symmetric,  (i.e.,  including  the  component  for 
(-a, -v)  if  the  coefficient  for  (u,v)  is  included,  we  can  be  sure  that  if /is  a  real-valued  image, 
then  each  of  the  images  in  the  decomposition  will  be  real-valued. 

Rotem  and  Zeevi  reconstruct  each  of  the  /^  from  the  information  sgn(/i),  k=],  ■  ■  ■  ,19. 
Since  each/^  is  a  bandpass  of/,  they  correspond,  crudely,  to  the  levels  of  the  scale-space  function 
V  that  we  use  to  reconstruct  /.  However,  the  analogy  is  poor,  since  the  /*  are  true  bandpass 
images,  with  one  octave  or  less  in  frequencies  in  at  least  one  of  the  two  spectral  dimensions, 
whereas  a  level  of  v  is  only  an  approximate  isotropic  bandpass  filter  of/.  Further,  the  number  of 
free  variables  in  any  given  /^  is  small  —  the  greatest  number  of  spectral  components  used  in  any 
of  the  fk  is  1024,  whereas  the  sgn(/t)  have  64K  bits  of  information,  for  each  k.  In  any  case,  the 
reconstruction  of  an/^  works  as  follows. 

Suppose  that  fk  is  bandpass  in  the  horizontal  spectral  dimension.  It  is  then  easily  shown 
that  each  row  of /*  is  a  one-dimensional  bandpass  signal.  (Alternatively,  if/*  is  bandpass  in  the 
vertical  spectral  dimension,  then  columns  of/t  are  also  bandpass).  The  spectral  support  of  the 
rows  or  columns  is  the  corresponding  band  in  the  horizontal  or  vertical  dimension  of  the  spectral 
domain  of/^.  The  reconstruction  begins  by  finding  individual  rows  or  columns  that  can  be  recon- 
structed from  the  sign  data  along  the  row  or  column.  For  a  one-dimensional  row  S  which  is 
bandpass,  denote  by  Bp  the  bandpass  filter,  thus  S=Bp(S).  Given  the  data  sgn(S),  let 
5o  =  Bp(sgn(5)).  The  algorithm  to  recover  S  proceeds  ileratively,  beginning  with  Sq,  and  itera- 
tively  setting 


>n+l 


=  S„ 


fi„(sgn(5„))-So 


The  iteration  is  not  guaranteed  to  converge,  but  if  it  does,  one  generally  has  that 
sgn(S„)  =  sgn(5).  Providing  all  the  conditions  for  uniqueness  of  the  representation  of  the  one- 
dimensional  signal  by  its  zero-crossing  are  met  (i.e.,  Logan's  conditions),  then  this  will  imply 
that  S„  and  S  differ  only  by  a  multiplicative  constant.  Having  reconstructed  a  number  of  rows 
(or  columns),  modulo  multiplicative  constants,  the  one-dimensional  signals  must  be  individually 
scaled.  By  reconstructing  a  one-dimensional  signal  in  a  transverse  direction,  ratios  can  be 
obtained  to  scale  the  individual  rows  (or  columns)  modulo  a  single  multiplicative  constant,  rela- 
tive to  the  whole  image.  If  the  image  is  bandpass  in  the  vertical  dimension,  then  columns  can  be 
reconstructed.  If  the  image  is  only  low-pass  in  this  dimension,  Zeevi  and  Rotem  show  that  it  will 
be  bandpass  in  a  diagonal  direction,  and  so  some  diagonal  can  be  reconstructed.  In  practice,  they 
recoastruct  several  transverse  one-dimensional  signals,  and  use  averages  to  best  scale  the  indivi- 
dual rows.  Reconstruction  along  columns  that  do  not  initially  converge  can  be  assisted  by  initiali- 
zation with  data  obtained  from  diagonals  or  rows,  and  from  horizontally  blurring  the  partially 
reconstructed  image. 

Once  each  /^  is  reconstructed  to  within  a  multiplicative  scale  factor,  a  scale  factor  is 
applied,  based  on  the  assumption  that  the  variance  of  each  f^  is  known  in  advance.  The  final 
reconstruction  of  /is  performed  using  the  reconstructed  ft,  the  stored  /o,  and  by  omitting  g.  The 
result  shown  in  their  paper  is  visually  good,  and  they  show  thai  errors  are  mostly  concentrated  on 
the  edges  of  the  original  image. 

3.3.  Sanz-Huang 

More  recently,  Jorge  Sanz.  and  Tom  Huang  have  shown  examples  of  reconstructions  of 
images  using  zero-crossings  of  the  scale-space  function  recorded  at  four  levels  [25]. 

First,  four  levels  of  the  scale-space  function  arc  chosen,  t,,  /=1,2,3,4,  where  v  is  the  .scale- 
space  function  of  the  initial  image  /    The  information  that  is  stored  in  the  representation  is 
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sgn(v(x,>',r,))  for  /=1,  •  ■  •  ,4.  The  operator  that  takes  f(x,y)  to  '£v(x,y,ti)  is  linear,  and  a  pseu- 
doinverse  is  computed.  This  operator,  which  cannot  restore  the  very  low  spectral  components  or 
the  very  high  spectral  components,  is  nonlheless  a  convolution  operator  (assuming  an  unbounded 
domain),  given  by  convolution  against  a  function  ^.  That  is, 

/  =  ({)*  IvC.-,r,). 

The  reconstruction  algorithm  proceeds  iteratively,  begirming  with 

So=(l)*Isgn(v(-,-.r,)) 

i 

Given  S„,  the  next  iteration  is  obtained  by  first  computing  the  scale-space  levels  of  the  estimate 

Sn. 

v„C,-,r,)  =  A^r(-,-,r,)*5„, 
and  then  modifying  v„  so  that  the  data  has  the  right  sign: 

J    \vn(x,y,t,)\   if  sgn(y(x,y,ti))>0 
\-\v„(x,y,ti)\   if  sgn(v(;c,y,f,))<0    " 

Finally,  the  v's  are  summed  and  filtered  to  obtain  the  next  estimate: 

5„+i=4>*Iv„(-,-.r,). 

This  procedure  is  not  guaranteed  to  converge,  but  was  used  for  experimental  purposes  any- 
way. When  applied  to  large,  natural  grayscale  images,  they  resulted  in  reconstructions  that  had 
scale-space  functions  whose  signs  agreed,  on  the  four  levels,  with  the  original  sign  data  of  the 
representation,  at  roughly  80%  of  the  pixels.  However,  the  reconstructions,  while  clearly  related 
to  the  originals,  had  a  rather  stylized  or  "impressionistic"  appearance,  and  many  perceptual 
differences  from  the  originals. 

3.4.  Comments 

Each  of  the  above  reconstruction  techniques  attempts  to  find  a  function  that  satisfies  the 
given  zero-crossing  information  by  taking  some  measurement  of  the  error  of  the  zero-crossings  of 
the  reconstruction  with  the  desired  zero-crossings.  We  might  call  this  generic  approach  one  of 
minimizing  data  error:  the  error  in  realizing  the  zero-crossing  data.  As  yet  another  version  of 
minimizing  data  error,  we  could  take  as  a  measure  of  the  error  the  square  integral  of  the  L^  differ- 
ence between  the  sign  of  the  scale-space  function  and  the  desired  sign  of  the  scale-space  function: 


IJ  II  sgn(v(-,-,0)-sgn(vo(-,-,0)  II  ^dt, 


where  sgn(vo)  is  the  given  data  in  scale-space,  and  v  =  AK*f.  This  functional  is  not  differentiable 
in  the  reconstructed  image  /due  to  the  discontinuous  signum  function.  However,  if  we  replace 
the  signum  function  sgn  with  a  C°°  approximate,  say  \)/,  then  the  error  measure  becomes  a  smooth 
function  of  the  /data.  Suppose  that  instead  of  being  given  the  sign  of  the  scale-space  data,  we  are 
instead  given  So(x,y,t)  =  '^(voix,y,t)).  Suppose  that  /is  the  current  estimate  of  the  original 
image.  Then  the  gradient  of  the  error  measure  with  respect  to  the  estimate  /is  proportional  to 


1 


[v/(v(-,-,f))-5o(-,-.r)]y(v(-,-,r))*A^: 


(x,y,t)  dt 


where  v  =AK*f.  Thus  a  simple  gradient  descent  procedure  can  be  formulated  to  minimize  the 
error.  The  result  is  a  particulariy  simple  linear  iterative  process.  Of  course,  the  representation  is 
now  the  approximate  signum  of  the  vq  data,  vj/(vo),  rather  than  the  simple  boolean  information 
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about  the  sign  of  the  filtered  data,  or  the  location  of  the  zero-crossings.  Further,  convergence  is 
not  guaranteed,  and  can  be  very  slow.  If  the  true  signum  data  is  used  as  the  target  data  in  place  of 
5o,  then  convergence  of  the  iterative  process  is  much  worse. 

Instead  of  pursuing  this  approach,  we  instead  develop  the  idea  of  minimizing  equation 
error,  discussed  in  the  next  section. 

4.  Formulation  and  Methods 

4.1.  Equation  error  minimization 

In  this  section,  we  describe  a  method  for  reconstructing  an  image  from  its  zero-crossings 
based  on  minimizing  equation  error.  The  method  was  reported  for  a  different  application  (deblur- 
ring)  in  conference  proceedings  [31],  and  preliminary  results  on  reconstructions  from  zero- 
crossings  were  shown  at  a  different  conference  [32]. 

Our  approach  to  reconstruction  from  zero-crossings  is  to  reconstruct  the  function  v  in 
scale-space,  and  make  use  of  the  fact  that  v  satisfies  the  Heat  Equation.  The  discrctized  data  can 
be  viewed  as  a  multilevel  grid  of  units  communicating  locally.  The  essence  of  the  idea  is  that  the 
units  should  achieve  values  satisfying  the  given  zero-crossing  constraints,  and  also  satisfying,  to 
the  extent  possible,  a  discrete  version  of  the  Heat  Equation.  As  in  a  network,  the  values  are 
updated  iteratively  using  information  from  local  values  to  minimize  a  measure  of  error. 

We  first  sketch  the  formulation  of  the  approach  in  a  continuous  domain.  Consider  the  ini- 
tial image  f{x)  and  its  scale-space  Laplacian-of-Gaussian  filtered  data  vo(;c,f).  We  assume  that 
the  zero-crossings  in  scale  space  of  vq  are  given,  in  the  form  of  knowledge  of  sgn[vo(x,f)]  at  all 
points.  We  then  pose  the  problem: 

Find  V  minimizing  11  Av  -  t:—  11  ^, 

dt 

subject  to  sgn[v(;c,r)]  =  sgn[vo(A:,r)]. 

Here,  the  norm  ||  •  ||  is  the  L^  norm  over  the  scale-space  R''xR'^. 

We  in  fact  advocate  a  slightly  different  measure  of  the  equation  error.  We  transform  the 
Heat  Equation,  a  second  order  partial  differential  equation,  into  a  first  order  system: 

Vv  =  o. 

Va=—  . 
dt 

The  problem  is  now  reformulated  as  the  search  for  a  pair  of  functions  (v,o)  minimizing 

l|V-o-|^||2-H||Vv-o||2 

subject  to 

sgn[v(A:.r)]  =  sgn[voU,r)]. 

4.2.  Discrete  formulation 

We  now  describe  the  discrete  formulation  of  these  ideas.  The  case  of  a  finite -difference 
discretization  for  n=\  will  be  described  in  some  detail.  Extension  to  two  spatial  dimensions 
involves  some  complications  which  will  be  summarized  more  briefly.  A  finite -clement  formula- 
tion is  also  given  for  the  ca.se  n  =  \. 
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From  the  given  initial  data/,,  we  recursively  construct  the  scale-space  function  v,  ^^  from  the 
formulas: 

4 

V,.i+1  =  —[v,-l,k+2Vi^k+Vi+l,kl 

4 

To  make  the  computation  feasible,  the  space  and  scale  domains  must  be  bounded.  Various 
approaches  to  handling  the  spatial  border  points  are  possible.  A  particularly  convenient  one, 
analogous  to  Neumann  boundary  conditions,  is  defined  as  follows.  Let  us  assume  that  /  is 
defined  for  0<i<N-\.  We  tlien  extend  /  to  be  periodic  with  period  2^,  according  to  the  formu- 
las: 

//v+,  =//v-]-,,     i=0,\,---  ,N-\, 

f2N+i=fn      all/. 

Then  the  values  Vi^  are  defined  for  all  /  and  all  k>0.  The  scale-space  function  v,^  is  periodic  with 
period  2N  in  the  variable  /,  and  has  the  same  structure  as  the  extended  /  data.  Specifically,  the 
data  from  i=0,  •  •  •  ,N-\  is  repeated  in  reverse  order  from  i=N,  ■  •  •  ,2N-\,  and  then  extended  as 
a  IN  block  periodically  (see  Figure  2  ). 

The  V  data  approximates  the  Laplacian  applied  to  the  Gaussian-blurred  image  data.  Associ- 
ated with  the  computed  scale-space  data  v,  ,t  is  the  sign  data  5,jfc,  which  is  defined  to  be  1, 0,  or  -1 
according  to  whether  v,jt  is  positive,  zero,  or  negative,  resf)ectively.  In  practice,  the  v  data  and 
resulting  5  data  are  calculated  only  for  0</  <A'-1,  and  for  Q<k<T,  for  some  arbitrary  cutoff  T. 

Having  computed  the  representation  in  the  form  of  the  S.k  data,  we  proceed  to  the  recon- 
struction. The  discrete  reformulation  of  the  Heat  Equation  as  a  first-order  system  involves  the 
introduction  of  a  second  grid  of  data,  <3,k  which  ideally  satisfies 

Then  the  equation  error  is  defined  by  the  quadratic  functional: 

v,-,t-v,-i,t 
::; O',* 


N-\ 
t>Oi=0 


v,,t+i-v,-A- 


<Si+\,k-^i.k 


We  minimize  E  subject  to  the  constraints 


->■  •*- 


T 


Figure  2.  Borders  are  handled  by  extending  the  original  data/.  The  data /o,  •  •  •  Jn-\  is  repeat- 
ed, alternately  from  left  to  right  and  then  from  right  to  left,  to  form  a  doubly  infinite  signal  that  is 
periodic  with  period  2A'.  Note  that  the  values  at  the  end  points,  e.g.,  the  datum  /ft-_, ,  appear  twice 
in  succession. 
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v,_i>0  when  S,,*  =  1 

Vii<0  whenS.jt  =  -1. 

A  further  constraint  is  needed  to  ensure  that  the  solution  v  =0  is  not  obtained.  In  principle,  this 
could  be  done  by  specifying  a  single  value  for  v  at  some  node.  In  practice  such  a  constraint 
would  propagate  only  very  slowly  through  the  grid,  making  convergence  difficult.  We  chose 
instead  to  specify  v  at  the  maximum  level  k=T,  where  it  is  ordinarily  quite  smooth.  Thus  little  is 
added  to  the  representation  by  inclusion  of  the  v,  7-  data.  (In  experiments  involving  fixing  gra- 
dient data  along  the  zero-crossings,  the  zero  solution  is  no  longer  admissible,  and  the  v,  7-  data  are 
not  included  in  the  representation.) 

In  the  experiments  described  below,  we  used  a  conjugate-gradient  approach  to  minimizing 
the  error  functional.  The  inequality  constraints  were  imposed  by  means  of  a  penalization  method, 
adding  to  £  a  term  that  grows  quadratically  in  the  variable  whenever  a  constraint  is  violated.  By 
placing  a  large  weight  on  the  penalization  term,  a  functional  suitable  for  minimizing  is  obtained 
that  will  ensure,  upon  convergence,  that  the  constraints  generally  hold.  Specifically,  we  use  a 
penalization  term  of  the  form: 

X~5:isgn(v,.*)-S,tl-v?t. 

The  gradient  of  the  equation  error  with  respect  to  the  unknowns  can  be  represented  as  a  sum 
of  discrete  convolutions  applied  to  the  data  v  and  o,  which  may  be  regarded  as  image  arrays. 
Then  for  E  equal  to  the  equation  error  as  above,  the  values  of  dE/dv,  ^  and  dE/do,  k  form  image 
arrays  also,  which  we  will  denote  by  V^E  and  V(,£.  The  formulas  (valid  for  A:>0)  are: 


Vv£  = 


Vo£  = 


0^0' 

'0   0   0" 

-1    10  -1 

*  V  -(- 

0^4 

0-4     0. 

.0     2  -2 

'-2     2  0' 

000' 

4^0 

*  V  -(- 

-1   6  -1 

0     0  0_ 

.000. 

*  a 


*  o 


Here,  the  three-by-three  arrays  are  masks  that  are  employed  in  three-by-three  local  convolutions 
against  the  array  data  v  and  c.  The  boundaries  can  be  conveniently  handled  by  placing  borders 
one  pixel  wide  around  the  data  arrays.  These  border  pixels  are  set  so  as  to  satisfy  the  chosen 
boundary  conditions.  The  sign  penally  adds  a  term 


A.lsgn(v)-5I 


•  V 


to  V„£. 


The  minimization  by  conjugate  gradients  [33]  is  implemented  as  follows.  Observe  that  the 
equation  error  E  as  fonmulated  above,  including  the  penally  terms  for  inequality  constraints,  is  a 
quadratic  functional  of  the  unknowns  v  and  o.  Combining  the  unknowns  into  a  single  vector,  say 
X,  the  equation  error  is  thus  equivalent  to 

E  =  ^x'^Ax-h'^x  +  Eo 

for  some  positive  definite  matrix  A,  vector  h.  and  scalar  constant  Eq.  Here  x^  denotes  transpose. 
Nonzero  values  for  llic  quantities  h  and  £0  arise  from  fixed  grid  values  on  the  boundaries.  The 
method  of  conjugate  gradients  is  ordinarily  formulated  as  the  minimization  of  a  functional  of 
exactly  this  form.  The  computation  involves  fonming  the  inner  products  of  certain  vectors  with 
each  other,  and  products  of  vectors  with  the  matrix  A.  However,  there  is  no  need  to  find  A,  h,  and 
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Eq  explicitly.  Instead,  we  will  make  use  of  the  fact  that  their  products  with  relevant  vectors  have 
known  forms.  Observe  that 

VE=Ax-h. 

Hence  if/?  is  an  arbitrary  vector  (e.g.  one  of  the  conjugate  gradients),  the  product  Ap  can  be  com- 
puted by  evaluating  the  expression  for  V£,  according  to  the  convolutions  given  above,  with  p 
substituted  for  x,  which  is  the  combined  v  and  o  data,  and  h  set  to  zero.  Setting  h  to  zero  is 
achieved  simply  by  replacing  all  constants  on  the  boundaries  by  zeroes.  This  also  has  the  effect 
of  setting  Eq  to  zero.  The  other  matrix  product  of  interest,  p^ Ap,  is  computed  by  evaluating 
directly  the  expression  for  £,  again  substituting  p  for  x  and  setting  boundary  constants  to  zero. 
Note  that  if  reflected  boundary  conditions  are  in  effect,  then  values  forp  must  be  reflected  in  the 
same  way.  Thus  we  see  that  all  the  matrix  products  necessary  for  a  conjugate  gradient  calcula- 
tion can  be  computed  without  knowing  A  explicitly. 

It  is  worth  noting  that  this  formulation  of  the  conjugate  gradient  minimization  is  readily 
adapted  to  heterogeneous  constraints  according  to  which  v  or  a  are  held  fixed  at  arbitrary  loca- 
tions in  the  grid.  One  simply  sets  the  values  of  p  to  0  at  the  corresponding  positions.  In  the  itera- 
tive updating,  those  points  are  skipped  over,  although  they  enter  into  the  calculation  of  the  equa- 
tion error  and  the  matrix  products. 

We  now  turn  to  the  finite-difference  discretization  for  the  case  «=2.  We  define  the  scale- 
space  grid  as  v,j  t,  ^>0,  and  formulate  the  Heat  Equation  as: 

Reformulation  of  this  equation  as  a  first-order  system  is  not  so  straightforward  as  in  the  one- 
dimensional  case.  A  fairiy  natural  approach  proceeds  with  introduction  of  variables  o{_y,t,  oj^^^jt, 
and  Wij^k  which  ideally  satisfy 

(1)  V..y.*-V/-lJ.<: 

g(2)^_  v..;.*-^.-.;-i.* 


and 

^i.j.k+\  ~  ^i.j.k  = 


^t+\,j,k       ^i,j,k 


^t,j+l,k      '^tj,k 


When  these  equations  hold  exactly,  the  data  v  is  blurred  from  level  k  to  level  A:-f-l  by  convolution 
against  the  two-dimensional  kernel  given  previously.  Analogously  to  the  n=\  case,  the  equations 
above  can  be  converted  to  a  quadratic  measure  of  equation  error  in  the  unknown  variables  v,  w, 
o^'\  and  o^^\  The  sign  constraints  on  v  can  again  be  implemented  using  quadratic  penalty  terms. 
The  resulting  error  functional  can  then  be  minimized  by  a  conjugate  gradient  procedure  in  the 
same  way  as  described  for  n-\.  It  is  worth  noting  that  this  approach  results  in  equations  for  the 
gradient  of  the  equation  error  which  are  linear  in  the  unknowns,  and  that  a  component  of  the 
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gradient  with  respect  to  an  unknown  at  (i,j,k)  depends  only  on  values  at  neighboring  grid  points. 
Thus  this  fonmulation  is  well  suited  to  a  hardware  implementation  involving  only  small-kernel 
convolutions  and  (for  the  conjugate  gradient  method)  accumulation  of  simple  global  sums. 

The  discretization  of  the  Heat  Equation  can  alternatively  proceed  by  a  fmite-element 
approach  [34].  For  the  case  n=\,  we  assume  a  regular  rectangular  discrete  grid,  iih,k-s), 
j=0,  •  ■  ■  N-\  and  k=0,\,  •  ■  •  .  Here  h  and  s  specify  the  scaling  of  the  two  axes  relative  to  each 
other.  (In  the  finite-difference  method  we  implicitly  set  /i=l/2,  5=1.)  The  simplest  elements  can 
be  formed  by  splitting  each  grid  square  into  two  triangular  regions,  along  one  of  the  two  chosen 
diagonals  (standard  triangular  elements);  see  Figure  3.  The  unknown  functions,  v  and  a,  can  be 
represented  as  continuous  piecewise  linear  functions,  linear  within  each  triangular  region,  by  sim- 
ply specifying  the  value  of  the  functions  at  each  grid  point.  Accordingly,  the  unknowns  are  the 
values  of  v(ih,ks)  and  C5(ih,ks),  for  all  /  and  k,  and  will  be  denoted  by  v, jt  and  o,  *.  The  simple 
piecewise  linear  elemenLs  suffice  since  the  equation  was  converted  to  a  first  order  system.  It  is 
then  a  simple  (albeit  tedious)  matter  to  compute  the  equation  error  for  such  functions  in  terms  of 
the  unknowns.  Like  the  finite-difference  formulation,  it  turns  out  that  the  equation  error  is  a  qua- 
dratic function  in  the  unknowns,  and  that  there  is  locality  of  dependence. 

The  case  /i=l/2,  s-l  is  very  common,  and  so  we  will  give  the  gradient  equations  for  this 
case.  As  was  done  for  the  finite-difference  case,  the  result  can  be  represented  as  a  sum  of  discrete 
convolutions,  with  the  data  v  and  o  considered  as  image  arrays.  Once  again,  except  near  border 
pixels,  these  gradient  values  3£/3v,  ^  and  dE/do,  k  can  be  viewed  as  images  obtained  by  convo- 
lutions against  the  grid  arrays  v,jt  and  o, _*.  The  formulas  are 
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Figure  3.  Triangular,  picccwisc-lincar,  finite  elements. 
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Separate  equations  are  needed  for  k  =  0  (9£/3v,  o  and  dEIdo,_o)  and  on  the  borders.  Note  that  a 
certain  amount  of  asymmetry  results  simply  from  choosing  triangular  elements  based  on  one  of 
two  possible  diagonal  cuts  through  each  rectangular  grid.  The  sign  penalty  adds  the  same  term  to 
V^E  as  for  the  finite-difference  formulation. 

4.3.  Image  reconstruction  from  scale-space  data 

Once  the  approximate  scale-space  function  v  has  been  reconstructed  by  the  above  minimi- 
zation process,  there  remains  the  problem  of  obtaining  the  image /to  which  it  corresponds.  First, 
the  discrete  reconstructed  data  is  summed  over  the  scale  dimension  to  give 

T 

The  r  data  can  be  seen  to  be  the  discrete  approximation  to  f-ur+i ,  where  the  «  data  is  defined  in 
terms  of  convolution  of /against  a  discrete  Gaussian  kernel,  as  discussed  in  section  2.2.2.  The 
problem  of  computing  /thus  reduces  to  inverting  I-Gj+i,  where  Gj+x  is  the  operator  yielding 
T+\  levels  of  Gaussian  blur.  The  simplest  and  cheapest  method  is  simply  to  incorporate  the  uj+i 
data  into  the  representation,  so  that  the  an  approximation  to  the  original  /  can  be  obtained  by 
computing  r  and  adding  ut+\  .  Since  the  ut+\  data  is  quite  blurry,  including  it  does  not  add  much 
information  to  the  representation. 

However,  since  the  operator/  -Gr+i  is  well-behaved,  it  is  possible  to  reconstruct/,  modulo 
an  additive  constant  and  a  multiplicative  constant,  by  computing  the  pseudoinverse  of /-G7+] 
using  a  singular- value  decomposition  (SVD).  The  resulting  /  always  has  an  average  value  of 
zero;  therefore  the  average  value  of  the  original  image  is  added  to  /  to  complete  the  reconstruc- 
tion. Unfortunately,  computing  the  SVD  can  be  prohibitively  expensive  for  images  larger  than 
about  16  by  16. 

Providing  the  operator  /-G7+1  can  be  viewed  as  a  convolution,  the  singular  value  decom- 
position can  be  computed  efficientiy,  even  for  large  images,  using  tiie  discrete  Fourier  transform 
of  /-G7+1,  or  equivalently,  the  inversion  may  be  done  by  means  of  the  discrete  convolution 
theorem.  Depending  upon  the  method  of  handling  borders,  /  -Gt+\  may  or  may  not  be  a  convo- 
lution. If  a  Fourier  method  is  applied  when  the  borders  are  handled  in  such  a  way  that  the  opera- 
tor is  not  a  convolution,  then  even  though  border  effects  are  confined  to  a  neighborhood  of  a  few 
pixels  from  the  border,  the  inversion  yields  poor  results.  However,  some  important  border 
methods  can  be  viewed  as  convolutions.  In  particular,  for  the  Neumann  boundary  conditions 
model,  as  discussed  in  the  previous  section,  the  operator  is  equivalent  to  a  convolution  operation 
against  a  doubly  periodic  image  that  is  twice  the  size  of  tiie  original  in  each  dimension,  obtained 
by  reflecting  the  image  in  each  dimension.  When  the  Fourier  method  is  used  to  invert  the  opera- 
tor/-Gr+i  computed  using  the  doubly-reflected  periodic  image,  Uie  results  are  excellent. 

The  inversion  is  complicated  by  tiie  singularity  that  arises  out  of  the  additive-constant 
ambiguity:  the  "DC-component"  of  I-Gt+\  has  zero  amplitude.  Further,  the  Fourier 
coefficients  belonging  to  low  frequencies  will  be  small  in  value.  A  cutoff  must  be  chosen  judi- 
ciously in  order  to  avoid  amplifying  low-frequency  errors  in  the  inversion  process. 
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Whether  the  operator  I  -Gj+i  is  a  convolution  or  not,  provided  it  is  linear,  and  provided  all 
eigenvalues  of  G  are  strictly  less  than  one,  the  inversion  may  be  obtained  from  the  series  expan- 
sion: 

(/-Gr+,r'=  ICGt-^iT. 

m=0 

which  can  usually  be  computed  quite  easily,  if  somewhat  intensively. 

The  results  obtained  by  using  the  appropriate  method  for  the  given  boundary  conditions  are, 
in  all  cases,  very  similar  to  those  obtained  using  the  simpler  method  of  adding  ut+\  to  r  as 
described  earlier,  assuming  that  ut+\  is  part  of  the  representation.  Consequently,  the  reconstruc- 
tions shown  below  were  performed  by  the  latter  method. 

5.  Results 

5.1.  One-dimensional  experiments 

Our  first  experiments  were  done  with  one-dimensional  data  taken  from  a  scanline  of  a  digi- 
tized image.  Figure  4  shows  a  plot  of  the  original  data.  The  Laplacian  monolith  vq  of  the  origi- 
nal image  was  computed  as  described  above,  to  a  depth  of  40  levels  of  blurring  by  the  "1  2  1" 
mask.  Figure  5  shows  the  sign  of  the  function  vq,  using  black  for  negative  regions,  white  for 
positive  regions,  and  gray  for  pixels  whose  values  lie  very  near  zero.  Reconstruction  proceeded 
by  minimization  of  equation  error,  using  conjugate  gradients  with  a  penalty  term  for  values  of  v 
which  differ  in  sign  from  the  corresponding  points  of  vq,  and  holding  the  topmost  (most  blurred) 
level  fixed  equal  to  the  corresponding  level  of  vq.  The  initial  estimate  for  v  was  sgn(vo),  and  the 
initial  estimate  foro  was  the  gradient  of  this  initial  v. 

The  computation  was  run  for  over  40,000  iterations,  achieving  an  equation  error  plus 
penalty  of  less  than  10"'°  with  only  77  sign  violations  out  of  5120  grid  points.    (It  should  be 


Figure  4.  A  scanline  taken  from  a  digitized  natural  image.  There  arc  128  pixels  in  the  scanline. 
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Figure  5.  The  sign  of  the  scale-space  filtered  data  shown  in  Figure  4.  Positive  regions  are 
shown  in  white,  negative  in  black,  and  values  very  near  zero  are  in  gray.  The  zero-crossing  oc- 
cur on  the  borders  between  black  and  white  regions,  and  along  contours  within  the  gray  regions. 


noted  that  the  sign  violations  are  small:  the  rms  average  value  of  the  violations  is  only  10"^  of  the 
rms  average  value  of  v.)  The  sum  of  levels  was  used  to  obtain  the  reconstructed  function,  by  the 
pseudoinverse  method  described  above.  This  reconstructed  function  is  shown  in  Figure  6,  and 
Figure  7  shows  its  sign  of  Laplacian  monolith. 

It  is  evident  that  the  zero-crossings  of  the  original  and  reconstructed  functions  are  essen- 
tially identical.  Also,  major  features  of  the  signal  are  reproduced,  although  there  are  significant 
differences  in  the  small-scale  details.    We  conclude  that  while  the  reconstruction  from  zero- 


Figure  6.  A  reconstruction  /  of  a  signal  from  the  reconstructed  Laplacian-of-Gaussian  filtered 
data  V.  This  function  should  be  the  same  as  Figure  4  if  stable  reconstruction  from  zero-crossings 
were  possible. 
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Figure  7.  The  sign  of  the  scale-space  filtered  data  shown  in  Figure  4.  The  data  is  nearly  identi- 
cal to  that  of  Figure  5,  showing  that  the  signals  in  Figures  4  and  6  have  virtually  the  same  zero- 
crossings  for  their  Laplacian-of-Gaussian  filtered  data  at  all  scales  of  resolution. 


crossings  preserves  some  of  the  essential  information  about  the  signal,  it  is  unstable  widi  respect 
to  high-frequency  components. 

Next,  the  reconstruction  process  was  applied  to  the  saine  data,  but  this  time  holding  fixed 
the  values  of  o  along  the  zero-crossings  of  vq,  and  leaving  free  the  values  of  v  at  maximum  blur. 
In  this  way,  we  have  specified  the  zero-crossings  and  the  gradient  information  at  the  zero- 
crossings  through  40  levels  in  scale  space.  After  approximately  6000  iterations,  the  equation 
error  plus  penalty  fell  below  10"'^.  The  reconstructed  function  is  not  shown,  because  it  is  visu- 
ally indistinguishable  from  the  original  signal.  The  reconstructed  function  differed  from  the  ori- 
ginal by  less  than  10~*  at  all  points  (on  an  intensity  scale  of  0  to  1).  This  result  can  be  taken  as 
empirical  evidence  that  the  representation  based  on  zero-crossings  plus  gradient  data  for  one- 
dimensional  images  is  complete  and  stable. 

5.2.  Two-dimensional  experiments 

We  next  turn  our  attention  to  two-dimensional  images.  The  original  images  are  64x64  pix- 
els with  8-bit  intensity  resolution.  Because  of  computing  constraints,  only  20  levels  of  scale  were 
used  in  the  reconstructions.  Figure  8  shows  an  original  image  of  a  pagoda,  and  the  result  of 
reconstruction  from  zero-crossings.  The  equation-error  minimization  was  run  for  approximately 
2000  iterations.  TTie  reconstruction  was  computed  by  subtracting  the  sum  of  levels  of  the  recon- 
structed Laplacian  monolith  from  the  blurred  original  image,  as  described  above.  As  can  be  seen, 
the  two  images  are  perceptually  different.  Mainly,  the  high-frequency  content  is  not  well  recon- 
structed, so  that  the  reconstruction  appears  as  a  slightly  blurred  version  of  the  original.  Despite 
the  evident  differences  between  the  two  images,  however,  the  zero-crossings  of  the  two  images 
are  nearly  identical.  There  are  sign  differences  at  less  than  2%  of  the  points  in  the  20-layer 
monolith,  concentrated  in  the  bonom  (high-frequency)  layers.  Vimjally  all  of  these  sign  differ- 
ences are  adjacent  to  zero-crossings,  causing  at  most  a  one-pixel  shift  in  the  location  of  the  zero- 
crossing. 

The  reconstruction  by  equation  error  minimization  is  sensitive  to  the  initial  guess  used  to 
start  the  conjugate-gradient  iteration.  For  the  reconstruction  in  Figure  8(b\  the  initial  Laplacian 
monolith  contained  values  of  ±1  everywhere,  with  the  sign  matching  the  required  sign  at  each 
point.  The  input  image  magnitude  was  in  the  range  [0,255].  Figure  9  shows  the  results  obtained 
using  initial  values  of  ±256  instead.  Even  after  runnmg  for  nearly  2(XXX)  iterations,  this  recon- 
struction shows  considerable  noise  and  artifacts.  TTiese  artifacts  take  the  form  of  strong  positive 
values  sidc-by-sidc  with  strong  negative  values,  which  largely  cancel  each  other  out  in  the 
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Figure  8.  An  original  image,  of  a  Japanese  pagoda,  and  a  reconstruction  from  the  zero-crossings 
in  scale-space. 
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Figure  9.  Another  reconstruction  of  the  pagoda  image  of  Figure  8(a),  using  a  different  initializa- 
tion. 


blurring  process.  Once  again,  the  zero-crossings  of  the  reconstructed  image  in  Figure  9  are  nearly 
the  same  as  those  of  the  original  image.   The  differences  between  the  reconstructed  images 
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(Figures  8(b)  and  9)  and  the  original  (Figure  8(a))  can  be  taken  as  an  empirical  demonstration  of 
the  lack  of  stability  of  the  representation  based  on  zero-crossings  alone. 

To  show  the  similarities  of  the  zero-crossing  representations,  we  show  the  sign  of  the 
Laplacian-of-Gaussian  for  the  original  image  and  the  reconstructed  images,  in  Figure  10,  at  a 
sampling  of  levels.  The  zero-crossings  arc  consequently  at  the  borders  between  black  and  white 
regions.  As  can  be  seen,  there  are  some  differences  between  the  zero-crossings,  especially  at  the 
bottom  levels.  However,  as  noted  before,  the  differences  consist  almost  exclusively  of  a  one- 
pixel  displacement  in  the  location  of  the  zero-crossings.  This  is  more  easily  seen  when  one 
displays  the  zero-crossings  as  curves,  and  flickers  between  the  two  images. 

It  is  worth  emphasizing  that  the  sign  of  Laplacian-of-Gaussian  in  Figure  10(b)  was  com- 
puted from  the  reconstructed  image,  and  is  not  the  sign  of  the  reconstructed  monolith.  The  solu- 
tion method  involves  the  reconstruction  of  a  function  in  scale-space,  v,  which  approximately 
satisfies  the  heat  equation,  subject  to  a  p)enalty  term  which  tends  to  make  sgn(v)  agree  with  the 
given  sgn(vo).  The  resulting  v,  summed  over  scale  space,  is  used  to  reconstruct  an  approximation 
/to  the  original  image/.  The  Laplacian-of-Gaussian  computed  from/,  call  it  vy,  should  approxi- 
mate vo,  the  Laplacian-of-Gaussian  of  the  original  image  /  However,  each  of  the  scale-space 
functions  vq,  v,  and  vy,  is  different.  It  is  the  sign  of  vy  which  is  shown  in  Figure  10,  representing 
the  true  locations  of  the  zero-crossings  of/  Accordingly,  Figures  8(a),  8(b),  and  9  represent 
three  images  whose  zero-crossings  are  nearly  the  same  at  all  levels,  and  yet  which  have  substan- 
tially different  numerical  and  perceptual  values.  Thus,  similarly  to  the  one -dimensional  case 
described  earlier,  these  results  empirically  establish  the  instability  of  the  zero-crossing  representa- 
tion, when  zero-crossings  are  used  alone.  Of  course,  one  should  modify  this  broad  statement  with 
the  note  that  a  particular  method  of  handling  borders  has  been  specified,  and  the  image  comes 
from  a  particular  class  of  images.  Also,  it  is  not  clear  how  much  the  reconstructions  might 
improve  if  the  iteration  were  allowed  to  run  for  much  longer  times. 

The  image  in  Figure  8(a)  contains  a  large  number  of  sharp  edges,  and  thus  high  frequency 
components.  It  is  known  that  in  the  conjugate-gradient  method  of  minimizing  equation  error,  the 
low-frequency  components  will  converge  before  the  high  frequencies.  In  Figure  11,  we  show  a 
standard  test  image.  The  image  has  significantly  smaller  Fourier  components  (nearly  a  factor  of 
10,  on  average)  at  the  higher  frequencies  as  compared  to  the  image  of  the  pagoda,  in  Figure  8. 
The  reconstruction  from  zero-crossings  of  this  image  works  quite  well,  as  shown  in  Figure  1 1(b). 
Again  there  are  in  fact  some  significant  numerical  differences  between  the  original  and  recon- 
struction in  Figure  11,  but  the  reconstructed  image  is  perceptually  similar  to  the  original,  except 
for  fine  details.  The  signs  of  the  Laplacian-of-Gaussians  of  the  images  arc  shown  in  Figurc  12. 
and  are  essentially  the  same,  as  in  Figure  10.  Once  again,  there  are  minor  differences  at  the  low 
levels  (high  frequencies). 

In  another  experiment,  not  shown,  using  an  even  softer  original  image  in  which  fine  detail 
was  relatively  lacking,  the  reconstructed  image  was  perceptually  almost  identical  to  the  original. 
Further,  the  iterative  process  for  reconstructing  v  converged  somewhat  more  quickly  to  a  low 
equation  error  with  this  image  than  with  the  image  of  the  pagoda. 

In  Figure  13,  we  show  the  result  of  a  reconstruction  of  Figure  8(a),  where  the  representation 
by  zero-crossings  has  been  enhanced  with  the  value  of  the  gradients  along  the  zero-crossing.  The 
reconstruction  is  now  stable,  accurate,  and  fast.  There  are  minor  numerical  differences,  but  no 
major  perceptual  differences  between  the  original  and  the  rcconstruction.  Of  course,  the  amount 
of  data  in  the  entire  rcprcseniation  exceeds  the  information  content  in  the  original  image,  since  at 
every  zero-crossing  point  on  every  level,  a  value  of  the  magnitude  of  the  gradient  of  the  filtered 
image  must  be  stored.  Since  these  values  are  stored  as  floating  point  numbers,  the  entire  data 
structure  is  in  no  way  a  compaction  of  the  original  image.  However,  a  compact  code  is  not  the 
objective  of  the  representation  by  zero-crossings,  so  that  at  this  juncturc,  we  can  conclude  that  a 
representation  that  includes  with  the  zero-crossings  the  gradient  data  along  the  zero-crossings  is 
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Figure  10.  Levels  0,  3,  6,  9,  12,  15,  and  18  of  the  sign  of  the  Laplacian  of  Gaussian  "monolith" 
of  Geft  column)  the  original  image,  (middle  column)  the  reconstructed  image  Figure  8(b),  and 
(right  column)  the  reconstructed  image  Figure  9. 
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Figure  11.  An  original  image,  and  the  reconstruction  using  the  zero-crossings.  The  original  im- 
age is  relatively  free  of  high  frequency  components. 
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Figure  12.  Levels  0,  3,  6,  9,  12,  15,  and  18  of  the  sign  of  the  Laplacian  monolith  of  (left  column) 
the  original  image,  and  (second  column)  the  reconstructed  image  of  Figure  1 1 . 


empirically  complete  and  stable. 

In  other  experiments  using  gradient  data,  we  have  found  that  reconstruction  is  stable  for  all 
kinds  of  images.  Noise  can  be  added  to  the  gradient  data,  and  stable  reconstruction  is  still  not 
difficult.  For  Figure  14,  multiplicative  noise  in  the  range  of  l±l/256  was  applied  to  the  gradient 
values  along  the  zero-crossings  before  reconstruction.  The  noise  simulates  storing  the  gradient 
values  with  8-bit  precision,  and  still  yields  a  good  representation.  However,  even  with  8-bit 
values  for  the  gradients,  the  total  amount  of  data  in  the  typical  representation  exceeds  the  data 
content  of  the  original  image. 

In  Figure  15,  we  show  the  same  image  reconstructed  from  a  representation  using  the  zero- 
crossings  together  with  10%  of  the  gradient  values.  The  locations  of  the  gradient  values  were 
chosen  at  random  from  the  set  of  all  zero-crossing  pixels.  Reconstruction  is  much  improved  from 
the  zero-crossings  alone.  However,  using  only  10%  of  the  gradient  values,  some  noise  is  present 
in  the  reconstruction.  If  the  percentage  of  gradient  values  retained  is  increased,  then  the  recon- 
struction is  improved.  Indeed,  with  50%  of  the  gradient  values,  reconstruction  is  approximately 
as  good  as  in  Figure  13,  where  100%  of  the  gradient  values  were  retained.  Good  reconstructions 
can  also  be  obtained  by  concentrating  the  gradient  data  in  the  higher- frequency  levels.  Figure  16 
shows  a  reconstruction  for  which  gradients  were  specified  on  all  zero-crossings  in  the  bonom  two 
levels,  and  nowhere  else.  The  number  of  zero-crossings  at  which  gradients  were  specified  came 
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Figure  13.  A  reconstruction  of  the  pagoda  image  (Figure  8(a))  using  zero-crossings  and  gradient 
data  along  the  zero-crossings. 
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Figure  14.  A  reconstruction  of  the  pagoda  image  (Figure  8(a))  using  zero-crossings  and  gradient 
data  along  zero-crossings,  with  the  gradients  contaminated  by  noise  of  relative  magnitude  1/256 
as  described  in  the  text. 
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Figure  15.  A  reconstruction  of  the  pagoda  image  (Figure  8(a))  using  zero-crossings  and  approxi- 
mately 10%  of  the  gradient  data  along  the  zero-crossings. 


to  about  22%  of  the  total.  The  reconstruction  is  nearly  as  good  as  that  for  which  gradients  were 
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completely  specified  at  all  levels.  This  result  shows  how  the  instability  is  greatest  for  the  high 
frequencies:  when  these  are  suitably  controlled,  the  reconstruction  proceeds  stably. 
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Appendix 

We  show  that  if  two  irreducible  polynomials  with  real  coefficients  in  n  real  variables  share 
in  common  the  same  (n-1 -dimensional)  zero-crossing  (even  for  just  an  open  patch),  then  the  two 
polynomials  are  the  same.  This  proof  was  supplied  by  Professor  D.  Mumford.  We  first  prove  an 
analog  of  the  decomposition  theorem  for  real  algebraic  varieties. 

Proposition  1:  Let  7:eR[xi,X2,  •  •  •  ,x„]  be  a  prime  ideal  and  let  /j  =  the  height  of  n,  (i.e.,  the 
length  of  a  maximal  chain  of  prime  ideals: 

(0)  =  Tto  ^  Tti  ^  7t2  5^  ■  ■  •  ^  rt^  =  7t. 

Then  if  V^(it)  is  the  locus  of  zeroes  in  R"  shared  by  all  polynomials  in  the  ideal  n,  (the  real 
variety  of  the  ideal  n)  then  Vjh(k)  is  the  union  of  real  submanifolds  of  R"  of  dimension  n-h  and 
less. 

Proof:  By  the  standard  decomposition  theorem,  the  complex  variety  V(n)  c  C,  which  is  the 
locus  of  complex  zeros  of  the  polynomials  in  n  viewed  as  polynomials  of  n  complex  variables, 
breaks  up  into  a  smooth  part,  which  is  always  an  «-/i  -dimensional  complex  submanifoid  of  C, 
and  a  singular  part,  which  is  a  finite  union  of  subsets  of  the  form  V(n),  with  n'  being  a  prime 
ideal  satisfying  Ti^rc',  (and  thus  n'  is  of  greater  height).  The  smooth  part  intersects  R"  in  an 
n-h  -dimensional  real  submanifoid  of  R",  and  so  the  proposition  follows  by  decreasing  induction 
on  the  height  of  the  prime  ideal. 

The  main  result  is: 

Proposition  2:  Let  peR[X],X2,  ' ' '  ,Xn]  be  an  irreducible  polynomial.  Thus  the  ideal  K  =  (p) 
generated  by  p  has  height  one.  Let  xe  V]r(p)  be  a  smooth  f)oint  (i.e.,  a  point  where  V^(p)  is  an 
n-\  dimensional  submanifoid),  and  let  i/  be  an  open  neighborhood  of  x  in  R".  TTien  for  any  real 
polynomial  q  of  n  variables,  if 

9=0  on  UnVjuin), 

then  p  divides  q. 

Proof:  Let  /  be  the  ideal  of  all  polynomials  in  ti.[x\,X2,  ■  ■  •  ,x„]  that  vanish  on  [/nVR(7i),  i.e., 
the  ideal  defined  by  this  subset.  The  radical  of  /,  denoted  y[T,  consists  of  all  polynomials  such 
that  some  power  of  the  polynomial  lies  in  /.  A  standard  theorem  in  algebraic  geometr>',  used  in 
the  proof  of  a  strong  form  of  the  Nullstellensatz,  shows  that  the  radical  of  /  is  the  finite  intersec- 
tion of  prime  ideals  containing  /,  so  that: 

V7"  =  Tt]  njt2  o  •  •  •  nn^,    ;:,  prime  ideals. 

Now,  the  patch  f/nViR(7t)  is  contained  in  Vf^{ylT),  by  the  definition  of/.  So, 

UnV^in)  c  Vr(n/7')  =  VR(7c,)uVR(n2)u  •  •  ■  uKR(rtt). 

Further,  pe/cVT,  so  that  peTt,  for  all  /.  Suppose  that  (p)g^  n,  for  evcr\'  /.  Then  the  height  of 
each  K,  is  at  lea.st  two,  and  thus  by  proposition  1,  the  variety  each  rR(7c,)  consists  of  submani- 
folds of  dimension  n-2  and  less.  Thus  the  union  of  the  VR(7t,)'s  cannot  contain  an  n-]  dimen- 
sional patch,  which  is  a  contradiction.  So  for  some  /.  we  must  have  (p)  =  n,,  and  .so 

(p)c/  cTi)  n  •  •■  nn*  ere,  =  (p). 
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Recalling  that  q  is  in  /,  we  have  therefore  that  q&  (p),  i.e.,  thatp  divides  q,  as  required.  ■ 

As  a  result  of  the  proposition,  we  have  a  uniqueness  theorem.  Suppose  that  p  and  q  are  two 
polynomials  with  domain  R".  Suppose  furthermore  that  Ap  and  Aq  are  irreducible  polynomials. 
The  corresponding  scale-space  functions  Vj  and  V2  will  also  be  polynomials,  (of  «+!  variables) 
and  must  be  irreducible,  for  otherwise  a  factorization  of  either  would  induce  (at  f=0)  a  factoriza- 
tion of  the  Laplacian  of  the  initial  polynomial.  Suppose  that  vi  and  V2  share  in  common  an  open 
patch  of  some  zero-crossing.  Recall  that  the  zero-crossings  are  defined  in  terms  of  the  boundaries 
of  the  regions  of  negativity  and  positivity,  and  thus  wiU  be  submanifolds  of  codimension  one. 
Applying  the  proposition  leads  to  the  conclusion  that  vj  and  V2  divide  one  another,  or 
equivalently,  are  equal.  Consequently,  p  and  q  are  equal. 

In  fact,  all  that  is  necessary  for  uniqueness  is  that  the  scale-space  functions  vj  and  V2 
should  be  irreducible.  To  extend  the  result,  it  would  be  interesting  to  characterize  all  factorable 
polynomials  that  satisfy  the  Heat  Equation. 
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